
% this file calculates welfare costs of losing Coal Job
% by iterating on value functions
%
% (1) Use estimation file for whole pop (Bestim) or cells (Bestimcells)
% (2) With (fut=0) or without (fut=1) allowing people to find jobs in coal sector 
%  

%% Real/Simulated data & transition/wage parameters %%
%
% from the data / from a simulation we need:
    Nx=length(param.lamNC);
    lamNC=param.lamNC;
    lamC=param.lamC;
    if fut==1
        disp('Calculating value of coal job, unemp & welfare costs of losing job with no coal jobs (post-coal-exit)')
        lamC=zeros(1,Nx).*lamC;
    elseif fut==0
        disp('Calculating value of coal job, unemp & welfare costs of losing coal job with coal jobs (pre-coal-exit).')
    end
    deltaNC=param.deltaNC;
    deltaC=param.deltaC;
    rho=param.rho;
    wCmu=param.wCmu;
    wCsig=param.wCsig;
    wNCmu=param.wNCmu;
    wNCsig=param.wNCsig;
%{
    wCmin = param.wCmin;    wCmax = param.wCmax;    wNCmin = param.wNCmin;    wNCmax = param.wNCmax;
%}
    %% numerical parameters
    Nx=length(wCmu); % for how many different groups are we calculating welfare losses?
    % when does algorithm stop?
    tol=0.0001;
    format shortG % avoid numerical answers in exponential format
    
    %% calibrated/set parameters
    % set discount rate - use monthly discount rate in what follows
    param.ry    = 0.1
    param.rm=param.ry/12;
    r=param.rm; %monthly
    param.bM=800;   % unemp benefits for periods 13+
                    % 800 euro for [ALG-2 + Wohngeld] (13+ months)
    % the intermediariy values bm(1,...,12) are set as function of wage
    % wCmedian=exp(wCmu);  % wCmean=exp(wCmu+sig2/2)

    
%% (1) initial guess for vector BM(x) 

% for each x, set initial BM guess: one third of median value of NC job
% (NB. exp(mu) is median of lognormal distribution)
BMinit=exp(wNCmu)/(3*r); %(monthly interest rate to discount monthly earnings)
% use guess to initialize algorithm
BMnu=BMinit;
BMold=ones(1,Nx);

%% (2) Calculate B_0 & then VC & VNC based on BM

algocount=0;
while abs(BMold-BMnu)>tol

% (2.1) Calculate B_0 based on updated (or initialized) BM
algocount=algocount+1;

BM=BMnu; % this is  valid independently of C or NC and xr...
% (intuition: when you're unemp long enough, your past work in C or NC 
%             no longer provides more than sweet memories, not extra cash)
%  (NB. this rests on our assumption that you don't retire out of unepm
%             otherwise different rho(xr) would create different BM(xr))
%
% Next step: Given long-run BM, calculate B0 depending on wage level and sector
% (xr & fut(=1 => no coal jobs to be found) determines value B0 by influencing likelihood of leaving
% unemp prior to period 13, where value BM hits)
B0=@(wr,xr)  Calgob(BM,param,wr,xr);

% (2.2) Calculate VNC & VC for different wage levels
%       (across distribution of w)
%
% Calgob is a function calculating value of becoming unemployed as f of past wage
VNCf=@(wr,xr) (1/(r + deltaNC(xr)+rho(xr))).*wr + deltaNC(xr).*B0(wr,xr);
VCf= @(wr,xr) (1/(r + deltaC(xr)+rho(xr))).*wr + deltaC(xr).*B0(wr,xr);

% min and max of observed mins/max over different Nx
%{
wCminmin=min(wCmin);   wCmaxmax=max(wCmax);
wNCminmin=min(wNCmin);wNCmaxmax=max(wNCmax);
%}
%% (3) Interpolate to get integrals in expression (9)

% can test whether expression (8) and  (9) give numerically the same?
VWNCF=@(wR,xr) VNCf(wR,xr).*lognpdf(wR,wNCmu(xr),wNCsig(xr));
VWCF =@(wR,xr) VCf(wR,xr).*lognpdf(wR,wCmu(xr),wCsig(xr));
% suploC=wCminmin; suphiC=wCmaxmax;  suploNC=wNCminmin; suphiNC=wNCmaxmax; suphiNC=exp(max(wNCmu).*(1+max(wNCsig)));
for nx=1:Nx
%{
    BMnu(nx)=1/(r+lamC(nx)+lamNC(nx)).*param.bM...
        +(lamNC(nx)/(r+lamC(nx)+lamNC(nx))).*integral(@(wR)VWNCF(wR,nx),suploNC,suphiNC)...
        +(lamC(nx)/(r+lamC(nx)+lamNC(nx))).*integral(@(wR)VWCF(wR,nx),suploC,suphiC);
%}
    BMnu(nx)=1/(r+lamC(nx)+lamNC(nx)).*param.bM...
        +(lamNC(nx)/(r+lamC(nx)+lamNC(nx))).*integral(@(wR)VWNCF(wR,nx),0,inf)...
        +(lamC(nx)/(r+lamC(nx)+lamNC(nx))).*integral(@(wR)VWCF(wR,nx),0,inf);
end
    fprintf('algorithm run number = %f\n', algocount)
    fprintf('new values of BM = %f\n', BMnu)
BMold=BM;
end

% Given model has been solved, can now estimate welfare costs deducting
% possibility of finding new job in Coal (B0fut)
% Now recalculate B0 without possibility of finding job in coal in future

% function handle for welfare costs of losing coal (equation (1)) 
WFC=@(wr,xr) VCf(wr,xr) - B0(wr,xr);


%% display results for different groups
matres=NaN(Nx,3);
makelabels=0;
if exist('rlab','var')==0
    makelabels=1;
    rlab=strings(size([1,Nx]));
end
clab = {'Value Coal Job', 'Value Unemp', 'Welfare Cost of Job Loss'};
for xn=1:Nx
    if makelabels==1 %check whether we have labelled the groups already in calling file
        disp('for group '); disp(xn)
        rlabx= num2str(xn);
        rlab(xn)= {['group' rlabx]}; % generic group names 
    end
    fprintf('mean value of coal job = %f\n', VCf(exp(wCmu(xn)+wCsig(xn)^2/2),xn))
    matres(xn,1) = VCf(exp(wCmu(xn)+wCsig(xn)^2/2),xn);
    
    fprintf('mean value of post-coal unemp = %f\n', B0(exp(wCmu(xn)+wCsig(xn)^2/2),xn))
    matres(xn,2) = B0(exp(wCmu(xn)+wCsig(xn)^2/2),xn);
    
    fprintf('mean welfare cost per person = %f\n',WFC(exp(wCmu(xn)+wCsig(xn)^2/2),xn)) 
    matres(xn,3) = WFC(exp(wCmu(xn)+wCsig(xn)^2/2),xn);
    res.wfcpp(xn)=WFC(exp(wCmu(xn)+wCsig(xn)^2/2),xn);
end

%% output results for latex
% adapt filename depending on which function calls this algorithm
path= '../../../paper/' ;
file=data.label;
filename= [path file];

% save results directly as latex table
matrix2latex(matres, filename, 'rowLabels', rlab, 'columnLabels', clab, 'alignment', 'c', 'format', '%-6.2f', 'size', 'small');

%{
fprintf('discounted prob of future bm streams from perspective of diff months= %f\n', discount)
fprintf('hypothetical discount prob without job-finding= %f\n', discountNOLAM)
fprintf('income differences moving from month 12 to 13= %f\n', incdiff(1))
%}

